{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Current and near-term quantum computers suffer from imperfections, as we repeatedly pointed it out. This is why we cannot run long algorithms, that is, deep circuits on them. A new breed of algorithms started to appear since 2013 that focus on getting an advantage from imperfect quantum computers. The basic idea is extremely simple: run a short sequence of gates where some gates are parametrized. Then read out the result, make adjustments to the parameters on a classical computer, and repeat the calculation with the new parameters on the quantum hardware. This way we create an iterative loop between the quantum and the classical processing units, creating classical-quantum hybrid algorithms.\n",
    "\n",
    "<img src=\"figures/hybrid_classical_quantum.svg\" alt=\"Hybrid classical-quantum paradigm\" style=\"width: 400px;\"/>\n",
    "\n",
    "These algorithms are also called variational to reflect the variational approach to changing the parameters. One of the most important example of this approach is the quantum approximate optimization algorithm, which is the subject of this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum approximate optimization algorithm\n",
    "\n",
    "The quantum approximate optimization algorithm (QAOA) is shallow-circuit variational algorithm for gate-model quantum computers that was inspired by quantum annealing. We discretize the adiabatic pathway in some $p$ steps, where $p$ influences precision. Each discrete time step $i$ has two parameters, $\\beta_i, \\gamma_i$. The classical variational algorithms does an optimization over these parameters based on the observed energy at the end of a run on the quantum hardware.\n",
    "\n",
    "More formally, we want to discretize the time-dependent $H(t)=(1-t)H_0 + tH_1$ under adiabatic conditions. We achieve this by Trotterizing the unitary. For instance, for time step $t_0$, we can split this unitary as $U(t_0) = U(H_0, \\beta_0)U(H_1, \\gamma_0)$. We can continue doing this for subsequent time steps, eventually splitting up the evolution to $p$ such chunks:\n",
    "\n",
    "$$\n",
    "U = U(H_0, \\beta_0)U(H_1, \\gamma_0)\\ldots U(H_0, \\beta_p)U(H_1, \\gamma_p).\n",
    "$$\n",
    "\n",
    "At the end of optimizing the parameters, this discretized evolution will approximate the adiabatic pathway:\n",
    "\n",
    "<img src=\"figures/qaoa_process.svg\" alt=\"Quantum approximate optimization algorithm\" style=\"width: 400px;\"/>\n",
    "\n",
    "The Hamiltonian $H_0$ is often referred to as the driving or mixing Hamiltonian, and $H_1$ as the cost Hamiltonian. The simplest mixing Hamiltonian is $H_0 = -\\sum_i \\sigma^X_i$, the same as the initial Hamiltonian in quantum annealing. By alternating between the two Hamiltonian, the mixing Hamiltonian drives the state towards and equal superposition, whereas the cost Hamiltonian tries to seek its own ground state.\n",
    "\n",
    "Let us import the necessary packages first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:09:25.405530Z",
     "start_time": "2018-11-19T20:09:25.393410Z"
    }
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "import numpy as np\n",
    "from functools import partial, reduce\n",
    "from qiskit import Aer, BasicAer, QuantumRegister, execute\n",
    "from qiskit.quantum_info import Pauli\n",
    "from qiskit.aqua import Operator\n",
    "from qiskit.aqua.components.initial_states import Custom\n",
    "from scipy.optimize import minimize\n",
    "np.set_printoptions(precision=3, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can define our mixing Hamiltonian on some qubits. As in the notebook on classical and quantum many-body physics, we had to define, for instance, an `IZ` operator to express $\\mathbb{I}\\otimes\\sigma_1^Z$, that is, the $\\sigma_1^Z$ operator acting only on qubit 1. We can achieve the same effect the following way (this time using the Pauli-X operator):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pauli_x(qubit, coeff):\n",
    "    eye = np.eye((n_qubits))\n",
    "    return Operator([[coeff, Pauli(np.zeros(n_qubits), eye[qubit])]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coefficient here means the strength of the transverse field at the given qubit. This operator will act trivially on all qubits, except the given one. Let's define the mixing Hamiltonian over two qubits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:08:17.863004Z",
     "start_time": "2018-11-19T20:08:17.843293Z"
    }
   },
   "outputs": [],
   "source": [
    "n_qubits = 2\n",
    "\n",
    "Hm = reduce(lambda x, y: x+y,\n",
    "            [pauli_x(i, 1) for i in range(n_qubits)])\n",
    "Hm.to_matrix()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, we will minimize the Ising problem defined by the cost Hamiltonian $H_c=-\\sigma^Z_1 \\otimes \\sigma^Z_2$. First let's create the functions defining the operators using the Pauli-Z matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pauli_z(qubit, coeff):\n",
    "    eye = np.eye((n_qubits))\n",
    "    return Operator([[coeff, Pauli(eye[qubit], np.zeros(n_qubits))]])\n",
    "\n",
    "\n",
    "def product_pauli_z(q1, q2, coeff):\n",
    "    eye = np.eye((n_qubits))\n",
    "    return Operator([[coeff, Pauli(eye[q1], np.zeros(n_qubits)) * Pauli(eye[q2], np.zeros(n_qubits))]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we define the cost Hamiltonian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:08:17.885938Z",
     "start_time": "2018-11-19T20:08:17.865184Z"
    }
   },
   "outputs": [],
   "source": [
    "J = np.array([[0,1],[0,0]])\n",
    "Hc = reduce(lambda x,y:x+y,\n",
    "        [product_pauli_z(i,j, -J[i,j])\n",
    "         for i,j in itertools.product(range(n_qubits), repeat=2)])\n",
    "Hc.to_matrix()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set $p=2$ and initialize the $\\beta_i$ and $\\gamma_i$ parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:08:17.925805Z",
     "start_time": "2018-11-19T20:08:17.905037Z"
    }
   },
   "outputs": [],
   "source": [
    "n_iter = 10 # number of iterations of the optimization procedure\n",
    "p = 2\n",
    "beta = np.random.uniform(0, np.pi*2, p)\n",
    "gamma = np.random.uniform(0, np.pi*2, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial state is a uniform superposition of all the states $|q_1,...,q_n\\rangle$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_state_vect = [1 for i in range(2**n_qubits)]\n",
    "init_state = Custom(n_qubits, state_vector=init_state_vect)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial circuit prepares the initial state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:08:18.058128Z",
     "start_time": "2018-11-19T20:08:18.050635Z"
    }
   },
   "outputs": [],
   "source": [
    "qr = QuantumRegister(n_qubits)\n",
    "circuit_init = init_state.construct_circuit('circuit', qr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a function `evolve` that takes a Hamiltonian $H$ and an angle $t$ and returns a circuit component made of the unitary matrix $e^{j H t}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:08:18.075525Z",
     "start_time": "2018-11-19T20:08:18.059933Z"
    }
   },
   "outputs": [],
   "source": [
    "def evolve(hamiltonian, angle, quantum_registers):\n",
    "    return hamiltonian.evolve(None, angle, 'circuit', 1,\n",
    "                              quantum_registers=quantum_registers,\n",
    "                              expansion_mode='suzuki',\n",
    "                              expansion_order=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create the circuit, we need to compose the different unitary matrice given by `evolve`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:08:18.091375Z",
     "start_time": "2018-11-19T20:08:18.077164Z"
    }
   },
   "outputs": [],
   "source": [
    "def create_circuit(qr, gamma, beta, p):\n",
    "    circuit_evolv = reduce(lambda x,y: x+y, [evolve(Hm, beta[i], qr) + evolve(Hc, gamma[i], qr)\n",
    "                                             for i in range(p)])\n",
    "    circuit = circuit_init + circuit_evolv\n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now create a function `evaluate_circuit` that takes a single vector `gamma_beta` (the concatenation of `gamma` and `beta`) and returns $\\langle H_c \\rangle = \\langle \\psi | H_c | \\psi \\rangle$ where $\\psi$ is defined by the circuit created with the function above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:08:18.117855Z",
     "start_time": "2018-11-19T20:08:18.093293Z"
    }
   },
   "outputs": [],
   "source": [
    "def evaluate_circuit(gamma_beta, qr, p):\n",
    "    n = len(gamma_beta)//2\n",
    "    circuit = create_circuit(qr, gamma_beta[:n], gamma_beta[n:], p)\n",
    "    return np.real(Hc.eval(\"matrix\", circuit, Aer.get_backend('statevector_simulator'))[0])\n",
    "evaluate = partial(evaluate_circuit, qr=qr, p=p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we optimize the angles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:08:37.699288Z",
     "start_time": "2018-11-19T20:08:18.119570Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: -0.9999999999999931\n",
       " hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>\n",
       "      jac: array([-0.,  0.,  0.,  0.])\n",
       "  message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'\n",
       "     nfev: 50\n",
       "      nit: 8\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([3.927, 6.113, 0.455, 2.749])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = minimize(evaluate, np.concatenate([gamma, beta]), method='L-BFGS-B')\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of the results\n",
    "\n",
    "We create a circuit using the optimal parameters found."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:08:37.795769Z",
     "start_time": "2018-11-19T20:08:37.701647Z"
    }
   },
   "outputs": [],
   "source": [
    "circuit = create_circuit(qr, result['x'][:p], result['x'][p:], p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the `statevector_simulator` backend in order to display the state created by the circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:08:38.175963Z",
     "start_time": "2018-11-19T20:08:37.797375Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.707 0.    0.    0.707]\n",
      "[-0.785  0.172  0.172 -0.785]\n"
     ]
    }
   ],
   "source": [
    "backend = BasicAer.get_backend('statevector_simulator')\n",
    "job = execute(circuit, backend)\n",
    "state = np.asarray(job.result().get_statevector(circuit))\n",
    "print(np.absolute(state))\n",
    "print(np.angle(state))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the state is approximately $e^{0.79j} \\frac{1}{\\sqrt{2}} \\left( |00 \\rangle + |11 \\rangle \\right)$. It corresponds to a uniform superposition of the two solutions of the classicial problem: $(\\sigma_1=1$, $\\sigma_2=1)$ and $(\\sigma_1=-1$, $\\sigma_2=-1)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now try to evaluate the operators $\\sigma^Z_1$ and $\\sigma^Z_2$ independently:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:08:38.181888Z",
     "start_time": "2018-11-19T20:08:38.178242Z"
    }
   },
   "outputs": [],
   "source": [
    "Z0 = pauli_z(0, 1)\n",
    "Z1 = pauli_z(1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:09:32.326007Z",
     "start_time": "2018-11-19T20:09:31.566128Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.7208456881689926e-15\n",
      "-1.7208456881689926e-15\n"
     ]
    }
   ],
   "source": [
    "print(Z0.eval(\"matrix\", circuit, Aer.get_backend('statevector_simulator'))[0])\n",
    "print(Z1.eval(\"matrix\", circuit, Aer.get_backend('statevector_simulator'))[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that both are approximatively equal to zero. It's expected given the state we found above and corresponds a typical quantum behavior where $\\mathbb{E}[\\sigma^Z_1 \\sigma^Z_2] \\neq \\mathbb{E}[\\sigma^Z_1] \\mathbb{E}[\\sigma^Z_2]$"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
